
@#define lags = (1:24)

var y_C_A a_A n_B_C n_B_D n_A_D n_A_C y_D_A  X_A C_A D_A L_A P_D_A P_C_A W_A pi_C_A pi_D_A g_A m y_C_B a_B y_D_B X_B C_B D_B L_B P_D_B P_C_B W_B pi_C_B pi_D_B g_B I_A I_B C_H_A C_H_B C_F_A C_F_B P_H_A P_H_B D_H_A D_F_A D_H_B D_F_B P_D_H_B P_D_H_A T B_A B_B T_A T_B P_A P_B Pi_A Pi_B i; 

varexo e;
   
parameters eta  psi alpha kappa kappa_N  beta delta rho mu_D_A mu_C_A mu_D_B mu_C_B phi omega phi_D tau;

// POSTERIOR MEAN, see Table 2 in paper
eta        =   1.0000000000;
phi        =   0.7000000000000;
omega      =   2.0000000000000;
psi        =   1.000000000000000;
beta       =   0.99500000000000;
alpha      =   0.25;
delta      =   0.05000000000;
mu_D_A     =   0.5000000000;
mu_C_A     =   0.50000000000;
mu_D_B     =   0.0;
mu_C_B     =   0.0;
phi_D      =   0.700;
rho        =   1.000000000000000;
tau        =   0.09 ;
kappa  =(1-alpha)*((1-beta*(1-delta))^(eta-1))/(((1-alpha)*((1-beta*(1-delta))^(eta-1)))+alpha) ;
kappa_N=((1-alpha)*(1-(1-delta)*beta)^(eta))/(delta*alpha+((1-alpha)*(1-(1-delta)*beta)^(eta))) ;

model(linear);
//production function

y_C_A = a_A + n_A_C ;
y_C_B = a_B + n_B_C ;

y_D_A = a_A + n_A_D ;
y_D_B = a_B + n_B_D ;

X_A = kappa*C_A + (1-kappa)*D_A ;
X_B = kappa*C_B + (1-kappa)*D_B ;

C_A = phi*C_H_A + (1-phi)*C_F_A ;
C_B = phi*C_H_B + (1-phi)*C_F_B ;

D_A = phi_D*D_H_A + (1-phi_D)*D_F_A ;
D_B = phi_D*D_H_B + (1-phi_D)*D_F_B ;

P_C_A = phi*P_H_A + (1-phi)*P_H_B ;
P_C_B = phi*P_H_B + (1-phi)*P_H_A ;

P_D_A = phi_D*P_D_H_A + (1-phi_D)*P_D_H_B ;
P_D_B = phi_D*P_D_H_B + (1-phi_D)*P_D_H_A ;

P_A = (1-alpha)*P_C_A + (alpha)*P_D_A;
P_B = (1-alpha)*P_C_B + (alpha)*P_D_B;

Pi_A = P_A - P_A(-1);
Pi_B = P_B - P_B(-1);

(((phi-1)/omega)- phi/eta)*C_H_A =(1-mu_C_A)*(((eta-1)/eta)*X_A - (1-phi)*((1/omega)-1/eta)*C_F_A +(L_A + P_H_A ))
@#for lag in lags 
   +(1-mu_C_A)*((mu_C_A)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_A - (1-phi)*((1/omega)-1/eta)*C_F_A +(L_A + P_H_A )) 
@#endfor
;


(((-phi)/omega)-(1-phi)/eta)*C_F_A =(1-mu_C_A)*(((eta-1)/eta)*X_A - ((1/omega)-1/eta)*phi*C_H_A + (L_A + P_H_B ))
@#for lag in lags 
   +(1-mu_C_A)*((mu_C_A)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_A - ((1/omega)-1/eta)*phi*C_H_A +  (L_A + P_H_B )) 
@#endfor
;


(((phi-1)/omega)- phi/eta)*C_H_B =(1-mu_C_B)*(((eta-1)/eta)*X_B - (1-phi)*((1/omega)-1/eta)*C_F_B + (L_B + P_H_B ))
@#for lag in lags 
   +(1-mu_C_B)*((mu_C_B)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_B - (1-phi)*((1/omega)-1/eta)*C_F_B + (L_B + P_H_B )) 
@#endfor
;



(((-phi)/omega)-(1-phi)/eta)*C_F_B =(1-mu_C_B)*(((eta-1)/eta)*X_B - ((1/omega)-1/eta)*phi*C_H_B + (L_B + P_H_A ))
@#for lag in lags 
   +(1-mu_C_B)*((mu_C_B)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_B -  ((1/omega)-1/eta)*phi*C_H_B + (L_B + P_H_A )) 
@#endfor
;






(((phi_D-1)/omega)-phi_D/eta)*D_H_A = (1-mu_D_A)*(((eta-1)/eta)*X_A + (1-phi_D)*((1/eta)-(1/omega))*D_F_A +(1/(1-beta*(1-delta)))*( P_D_H_A + L_A - beta*(1-delta)*(P_D_H_A(+1) +L_A(+1) ))) 
@#for lag in lags 
+ (1-mu_D_A)*((mu_D_A)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_A + (1-phi_D)*((1/eta)-(1/omega))*D_F_A +(1/(1-beta*(1-delta)))*( P_D_H_A + L_A - beta*(1-delta)*(P_D_H_A(+1) +L_A(+1) ))) 
@#endfor
;



(((-phi_D)/omega)-(1-phi_D)/eta)*D_F_A  = (1-mu_D_A)*(((eta-1)/eta)*X_A + (phi_D)*((1/eta)-(1/omega))*D_H_A +(1/(1-beta*(1-delta)))*( P_D_H_B + L_A- beta*(1-delta)*(P_D_H_B(+1) +L_A(+1) ))) 
@#for lag in lags 
+ (1-mu_D_A)*((mu_D_A)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_A + (phi_D)*((1/eta)-(1/omega))*D_H_A +(1/(1-beta*(1-delta)))*( P_D_H_B + L_A- beta*(1-delta)*(P_D_H_B(+1) +L_A(+1) ))) 
@#endfor
;



(((phi_D-1)/omega)-phi_D/eta)*D_H_B = (1-mu_D_B)*(((eta-1)/eta)*X_B + (1-phi_D)*((1/eta)-(1/omega))*D_F_B +(1/(1-beta*(1-delta)))*( P_D_H_B + L_B - beta*(1-delta)*(P_D_H_B(+1) +L_B(+1) ))) 
@#for lag in lags 
+ (1-mu_D_B)*((mu_D_B)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_B + (1-phi_D)*((1/eta)-(1/omega))*D_F_B +(1/(1-beta*(1-delta)))*( P_D_H_B + L_B - beta*(1-delta)*(P_D_H_B(+1) +L_B(+1) ))) 
@#endfor
;



(((-phi_D)/omega)-(1-phi_D)/eta)*D_F_B  = (1-mu_D_B)*(((eta-1)/eta)*X_B + (phi_D)*((1/eta)-(1/omega))*D_H_B +(1/(1-beta*(1-delta)))*( P_D_H_A + L_B- beta*(1-delta)*(P_D_H_A(+1) +L_B(+1) ))) 
@#for lag in lags 
+ (1-mu_D_B)*((mu_D_B)^(@{lag}))*EXPECTATION(-@{lag})(((eta-1)/eta)*X_B + (phi_D)*((1/eta)-(1/omega))*D_H_B +(1/(1-beta*(1-delta)))*( P_D_H_A + L_B- beta*(1-delta)*(P_D_H_A(+1) +L_B(+1) ))) 
@#endfor
;


//(1/eta)*((1-eta)*(X_A - X_A(-1)) - C_A + C_A(-1)) - (1/eta)*((1-eta)*(X_B - X_B(-1)) - C_B + C_B(-1))  = (pi_C_A - pi_C_B) ;
L_A=L_B ;

//(1/eta)*((1-eta)*(X_A ) - C_A ) - (1/eta)*((1-eta)*(X_B ) - C_B )  = (P_C_A - P_C_B) ;

psi*(kappa_N*n_A_C + (1-kappa_N)*n_A_D) = (-(tau/(1-tau))*T_A + W_A + L_A)  ;
psi*(kappa_N*n_B_C + (1-kappa_N)*n_B_D) = (-(tau/(1-tau))*T_B + W_B + L_B) ;


pi_C_A= P_C_A - P_C_A(-1) ;
pi_C_B= P_C_B - P_C_B(-1) ;

pi_D_A= P_D_A - P_D_A(-1) ;
pi_D_B= P_D_B - P_D_B(-1) ;


P_H_A  = W_A - a_A ;

P_H_B  = W_B - a_B ;

P_D_H_A  = W_A - a_A ;

P_D_H_B  = W_B - a_B ;


y_C_A = g_A + phi*C_H_A + (1-phi)*C_F_B ;
y_C_B = g_B + phi*C_H_B + (1-phi)*C_F_A ;


delta*I_A = (D_A-(1-delta)*D_A(-1)) ;
delta*I_B = (D_B-(1-delta)*D_B(-1)) ;

y_D_A = g_A + phi_D*(1/delta)*(D_H_A-(1-delta)*D_H_A(-1)) + (1-phi_D)*(1/delta)*(D_F_B-(1-delta)*D_F_B(-1)) ;
y_D_B = g_B + phi_D*(1/delta)*(D_H_B-(1-delta)*D_H_B(-1)) + (1-phi_D)*(1/delta)*(D_F_A-(1-delta)*D_F_A(-1)) ;

T_A=  T ;
T_B=  T ;

m = 0.5*(kappa_N*P_H_A + kappa_N*y_C_A + (1-kappa_N)*P_D_H_A + (1-kappa_N)*y_D_A ) + 0.5*(kappa_N*P_H_B + kappa_N*y_C_B + (1-kappa_N)*P_D_H_B + (1-kappa_N)*y_D_B ) ;

kappa_N*(P_C_A(+1) + C_A(+1)) + (1-kappa_N)*(P_D_A(+1) + I_A(+1)) + B_A(+1) = (W_A(+1) + (kappa_N*n_A_C(+1) + (1-kappa_N)*n_A_D(+1))) + (1/beta)*B_A + i  ; 
 
kappa_N*(P_C_B(+1) + C_B(+1)) + (1-kappa_N)*(P_D_B(+1) + I_B(+1)) + B_B(+1) = (W_B(+1) + (kappa_N*n_B_C(+1) + (1-kappa_N)*n_B_D(+1))) + (1/beta)*B_B + i ;

B_A = -B_B ;

T=rho*T(-1) +e(-8) ;
//T=0;
m = rho*m(-1)  ;

g_A = rho*g_A(-1);
g_B = rho*g_B(-1) ;

a_A=0;
a_B=0;
//a_A=rho*a_A(-1) +e(-8);
//a_B=rho*a_B(-1) +e(-8);

// shock processes

end;

//check;
steady;

shocks;
var e;      stderr 1.0;
//var xx;      stderr 0.06;
//var zz;      stderr 0.06;
end;



stoch_simul( order=1,irf=24, nograph);


//stoch_simul(order=1,irf=32);
//stoch_simul;
